##########################################################################################

library(data.table)
library(optparse)
library(ArchR)
library(ggsci)
library(tidyverse)
library(BSgenome.Hsapiens.UCSC.hg38)
if(1!=1){
  library(TFBSTools)
  library(JASPAR2020)
}

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--peak_gene_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata"

    ## peak-gene文件
    peak_gene_file <- "~/20231121_singleMuti/results/celltype_plot/peak2gene/germ/peakToGeneHeatmap_LabelClust_k25.tsv"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/qc_atac_v3/germ_peak-gene"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
peak_gene_file <- opt$peak_gene_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(comine_data_file)
## testis_combined_peak_combineRNA
dat_peak_gene <- data.frame(fread(peak_gene_file))

###########################################################################################
## 输出到对应目录
projHeme5 <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = rownames(testis_combined_peak_combineRNA@cellColData),
  outputDirectory = out_path,
  dropCells = TRUE,
  force = TRUE
)

###########################################################################################
## 提取感兴趣的peak集合
tmp_peak <- data.frame(projHeme5@peakSet)
tmp_peak <- paste0(tmp_peak$seqnames , ":" , tmp_peak$start , "-" , tmp_peak$end)
use_peak <- tmp_peak[tmp_peak %in% unique(dat_peak_gene$peak)]
tmp_peak_index <- which(tmp_peak %in% use_peak)

## peak只用存在peak-gene关系的peak
projHeme5@peakSet <- projHeme5@peakSet[tmp_peak_index,]
## 使用自己定义的peak集合
projHeme5 <- addPeakSet(projHeme5,projHeme5@peakSet,force = TRUE)
projHeme5 <- addPeakMatrix(projHeme5)

## 重新计算这部分peak上的motif活性，背景用peak-gene的peak

if(1!=1){
  #pwm_jaspar <- TFBSTools::getMatrixSet(
  #  x = JASPAR2020::JASPAR2020,
  #  opts = list(
  #  tax_group = "vertebrates",
  #  species = "Homo sapiens",
  #  matrixtype = "PWM"
  #  )
  #)

  #pwm_set <- TFBSTools::getMatrixSet(x = JASPAR2020::JASPAR2020, opts = list(all_versions = FALSE, species = 9606, collection="CORE", matrixtype="PWM"))
  #pwm_set2 <- TFBSTools::getMatrixSet(x = JASPAR2020::JASPAR2020, opts = list(all_versions = FALSE, species = 9606, collection="UNVALIDATED", matrixtype="PWM"))
  #pwm_set <- c(pwm_set1, pwm_set2)
  #names(pwm_set) <- name(pwm_set) #converts Motif IDs to gene names
  #pwm_set[duplicated(names(pwm_set)),] #identifies duplicates
  #names(pwm_set) <- make.unique(names(pwm_set)) #adds a unique suffix onto end of duplicates
  #names(pwm_set) <- sapply( strsplit(sapply(strsplit(names(pwm_set) , "[(]") , "[" , 1) , ":" ) , "[" , 1)

  #projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "JASPAR2020" , motifPWMs = pwm_jaspar , name = "Motif",force = TRUE)
}
projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "Cisbp" , name = "Motif",force = TRUE)
projHeme5 <- addBgdPeaks(projHeme5 , force = TRUE)
projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "Motif",
  force = TRUE
)

projHeme5 <- addPeak2GeneLinks(ArchRProj = projHeme5 , useMatrix = "GeneExpressionMatrix")


###########################################################################################
## 输出
testis_combined_peak_combineRNA <- projHeme5
image_name <- paste0( out_path , "/testis_combined_peak.combineRNA.qc.Rdata")
save(  testis_combined_peak_combineRNA , file = image_name )
saveArchRProject( testis_combined_peak_combineRNA , out_path )

